library(tidyverse)
library(plotly)
library(sf)
library(mapview)
library(tigris)
library(censusapi)
library(leaflet)
library(lehdr)


options(
  tigris_class = "sf",
  tigris_use_cache = TRUE
)

Sys.setenv(CENSUS_KEY="10dcd73d7c043e91bac9fb8d3989cbff54b08790")

Load social distancing data and blockgroups

Load the Safegraph social distancing data and San Jose blockgroups

# get SJ blockgroups 
# get San Jose block groups
scc_blockgroups <- block_groups("CA","Santa Clara", cb=F, progress_bar=F)

# Use tracts sent to us by San Jose
sj_tracts <- st_read("/Users/simonespeizer/pCloud Drive/Shared/SFBI/Data Library/San_Jose/CSJ_Census_Tracts/CSJ_Census_Tracts.shp") %>%
  st_as_sf() %>%
  st_transform(st_crs(scc_blockgroups))
## Reading layer `CSJ_Census_Tracts' from data source `/Users/simonespeizer/pCloud Drive/Shared/SFBI/Data Library/San_Jose/CSJ_Census_Tracts/CSJ_Census_Tracts.shp' using driver `ESRI Shapefile'
## Simple feature collection with 219 features and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 6112856 ymin: 1869687 xmax: 6255982 ymax: 1996555
## epsg (SRID):    2227
## proj4string:    +proj=lcc +lat_1=38.43333333333333 +lat_2=37.06666666666667 +lat_0=36.5 +lon_0=-120.5 +x_0=2000000.0001016 +y_0=500000.0001016001 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs
sj_citycouncil_disticts <- st_read("/Users/simonespeizer/pCloud Drive/Shared/SFBI/Data Library/San_Jose/City Council Districts/CITY_COUNCIL_DISTRICTS.shp") %>%
  st_as_sf() %>%
  st_transform(st_crs(scc_blockgroups))
## Reading layer `CITY_COUNCIL_DISTRICTS' from data source `/Users/simonespeizer/pCloud Drive/Shared/SFBI/Data Library/San_Jose/City Council Districts/CITY_COUNCIL_DISTRICTS.shp' using driver `ESRI Shapefile'
## Simple feature collection with 10 features and 7 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 6112856 ymin: 1869687 xmax: 6255982 ymax: 1996555
## epsg (SRID):    2227
## proj4string:    +proj=lcc +lat_1=38.43333333333333 +lat_2=37.06666666666667 +lat_0=36.5 +lon_0=-120.5 +x_0=2000000.0001016 +y_0=500000.0001016001 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs
# from code written by others to get SJ blockgroups
sj_blockgroups <-
  scc_blockgroups %>%
  st_centroid() %>%
  st_join(sj_tracts, left = F) %>%
  st_join(sj_citycouncil_disticts%>% dplyr::select(DISTRICTS)) %>%
  mutate(
    DISTRICTS = DISTRICTS %>% factor(levels = c("1","2","3","4","5","6","7","8","9","10"))
  ) %>%
  st_set_geometry(NULL) %>%
  left_join(scc_blockgroups%>% dplyr::select(GEOID), by = "GEOID") %>%
  st_as_sf() %>%
  dplyr::select(GEOID, DISTRICTS)

# the spatial join leaves off two blockgroups which are touching district 9. The following code assigns those to district 9
sj_blockgroups$DISTRICTS[is.na(sj_blockgroups$DISTRICTS)] <- 9

# code from others in the class to get social distancing data 
sj_socialdistancing <- readRDS("/Users/simonespeizer/pCloud Drive/Shared/SFBI/Restricted Data Library/Safegraph/covid19analysis/sj_socialdistancing.rds") %>% 
  mutate(date = date_range_start %>%  substr(1,10) %>% as.Date()) %>% 
  left_join(sj_blockgroups, by = c("origin_census_block_group" = "GEOID")) %>% 
  filter(!is.na(DISTRICTS))

# obtaining weekends
weekends <-
  sj_socialdistancing %>% 
  filter(!duplicated(date)) %>% 
  arrange(date) %>% 
  mutate(
    weekend = 
      ifelse(
        (date %>% as.numeric()) %% 7 %in% c(2,3),
        T,
        F
      )
  ) %>% 
  dplyr::select(date,weekend)

sj_socialdistancing <- 
  sj_socialdistancing %>% 
  left_join(weekends)

# date of the shelter in place order
shelter_start <- "2020-03-16" %>% as.Date()

# get average staying at home on weekdays in January and February
sj_pre_sd_at_home_average <- sj_socialdistancing %>% 
  filter(weekend == F) %>% 
  filter(date <  as.Date("2020-03-01")) %>%
  group_by(origin_census_block_group) %>% 
  summarize(completely_home_device_count = sum(completely_home_device_count), device_count = sum(device_count)) %>% 
  mutate(`% Completely at Home Pre Shelter` = (completely_home_device_count/device_count*100) %>% round(1), `% not completely at home pre shelter` = (100 - `% Completely at Home Pre Shelter`))

Obtaining demographic variables

Here I obtain various demographic data, including income (percent below 50% and 80% of area median income), vehicle ownership, age, English language ability, and occupants per room.

# obtain the saved census data 
setwd("~/Documents/2020 Spring Quarter/CEE 218Z")
acs_vars = readRDS("censusData2018_acs_acs5.rds")
setwd("~/Documents/2020 Spring Quarter/CEE 218Z/covid19")
# load in income data - code adapted from other students
sj_median_income_by_block <-
  getCensus(
    name = "acs/acs5",
    vintage = 2018,
    region = "block group:*", 
    regionin = "state:06+county:085",
    vars = "B19013_001E"
  ) %>%
  mutate(
    blockgroup =
      paste0(state,county,tract,block_group)
  ) %>% 
  select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>% 
  rename(
    Median_Income = B19013_001E 
  ) %>% 
  filter(!is.na(Median_Income)) %>% 
  left_join(sj_blockgroups, by = c("blockgroup" = "GEOID")) %>% #this code gives each blockgroup a district designation
  filter(
    !is.na(DISTRICTS)
  ) %>% 
  
  # this code joins our census data with the social distancing data, processed as shown below
  left_join(sj_socialdistancing %>%  
                          filter(weekend == F) %>% 
                          filter(date > shelter_start) %>%
                          group_by(origin_census_block_group) %>% 
                          summarize(
                                    completely_home_device_count = sum(completely_home_device_count),
                                    device_count = sum(device_count)) %>% 
                          mutate(`% Completely at Home` = (completely_home_device_count/device_count*100) %>% round(1), 
                                 `% not completely at home` = (100 - `% Completely at Home`)),
            by = c("blockgroup" = "origin_census_block_group")
  ) %>% 
  filter(
    !is.na(device_count)
  ) %>% 
  left_join(sj_pre_sd_at_home_average %>% dplyr::select(origin_census_block_group, `% Completely at Home Pre Shelter`, `% not completely at home pre shelter`), by = c("blockgroup" = "origin_census_block_group"))

sj_ami_by_block <-
  getCensus(
    name = "acs/acs5",
    vintage = 2018,
    region = "block group:*", 
    regionin = "state:06+county:085",
    vars = "group(B19001)"
  ) %>%
  mutate(
    blockgroup =
      paste0(state,county,tract,block_group)
  ) %>% 
  select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>% 
  dplyr::select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
  group_by(blockgroup) %>% 
  summarize(
    Total = B19001_001E,
    `Under 75,000` = sum(B19001_002E, B19001_003E, B19001_004E, B19001_005E, B19001_006E, B19001_007E, B19001_008E, B19001_009E, B19001_010E, B19001_011E, B19001_012E),
    #sum(lapply(2:12, function(x) as.name(paste0("B19001_00",x,"E"))))
    `Under 100,000` = sum(B19001_002E, B19001_003E, B19001_004E, B19001_005E, B19001_006E, B19001_007E, B19001_008E, B19001_009E, B19001_010E, B19001_011E, B19001_012E, B19001_013E), 
    `Under 125,000` = sum(B19001_002E, B19001_003E, B19001_004E, B19001_005E, B19001_006E, B19001_007E, B19001_008E, B19001_009E, B19001_010E, B19001_011E, B19001_012E, B19001_013E, B19001_014E)
  ) %>% 
  mutate(
    `% under 75,000` = `Under 75,000` / Total * 100,
    `% over 75,000` = (100 - `% under 75,000`),
    `% under 100,000` = `Under 100,000` / Total * 100,
    `% over 100,000` = (100 - `% under 100,000`),
    `% under 125,000` = `Under 125,000` / Total * 100,
    `% over 125,000` = (100 - `% under 125,000`),
  ) %>% 
  left_join(sj_median_income_by_block %>% dplyr::select(-Median_Income)
  ) %>% 
  filter(!is.na(device_count))
# loading in language data - code adapted from other students
sj_lang_by_block <-
  getCensus(
    name = "acs/acs5",
    vintage = 2018,
    region = "block group:*", 
    regionin = "state:06+county:085",
    vars = "group(B16004)"
  )  %>% 
  mutate(
    blockgroup =
      paste0(state,county,tract,block_group)
  ) %>% 
  select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>% 
  select(-c(contains("EA"),contains("MA"),contains("M"))) %>% 
  gather(
    key = "variable",
    value = "estimate", 
    - blockgroup
  ) %>% 
  left_join(acs_vars, by = c("variable" = "name")) %>% 
  mutate(
    tier = substr(label,lapply(label, function(x) max(unlist(gregexpr('!!',x)))+2),nchar(label))
  ) %>% 
  filter(tier %in% c('Speak English "not well"', 
                     'Speak English "not at all"', 
                     'Total', 'Speak Spanish', 
                     'Speak Asian and Pacific Island languages')) %>% 
  group_by(blockgroup, tier) %>% 
  summarise(
    estimate1 = sum(estimate)
  ) %>% 
  spread(
    key = "tier",
    value = "estimate1"
  ) %>% 
  mutate(
    `% speaking english < well` = (`Speak English "not well"` + `Speak English "not at all"`) / Total * 100,
    `% speaking english > well` = (100 - `% speaking english < well`),
    `% speaking spanish` = (`Speak Spanish`/ Total) * 100,
    `% not speaking spanish` = (100 - `% speaking spanish`),
    `% speaking api` = (`Speak Asian and Pacific Island languages` / Total) * 100
  ) %>% 
  left_join(sj_median_income_by_block %>% dplyr::select(-Median_Income)) %>% 
  filter(!is.na(device_count)) %>% 
  mutate(log_perc = log(`% speaking english < well`))
# loading in age data - specifically looking at percentage 65+ and percentage <30
sj_age_by_block <- getCensus(
    name = "acs/acs5",
    vintage = 2018,
    region = "block group:*", 
    regionin = "state:06+county:085",
    vars = "group(B01001)"
  ) %>% 
  mutate(
    blockgroup =
      paste0(state,county,tract,block_group)
  ) %>% 
  select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>% 
  select(-c(contains("EA"),contains("MA"),contains("M"))) %>% 
  gather(
    key = "variable",
    value = "estimate", 
    - blockgroup
  ) %>% 
  mutate(
    label = acs_vars$label[match(variable,acs_vars$name)]
  ) %>% 
  select(-variable) %>% 
  separate(
    label,
    into = c(NA,NA,"sex","age"),
    sep = "!!"
  ) %>% filter(!is.na(age)) %>% 
  mutate(elderly = ifelse(age %in% c("65 and 66 years", "67 to 69 years", "70 to 74 years", "75 to 79 years", "80 to 84 years", "85 years and over"), estimate, NA), `less than 30` = ifelse(age %in% c("Under 5 years", "5 to 9 years", "10 to 14 years", "15 to 17 years", "18 and 19 years", "20 years", "21 years", "22 to 24 years", "25 to 29 years"), estimate, NA)) %>% 
  group_by(blockgroup) %>% 
  summarize(elderly = sum(elderly, na.rm = T), `less than 30` = sum(`less than 30`, na.rm = T), total = sum(estimate, na.rm = T)) %>% 
  mutate(`percent elderly` = elderly*100 / total, `percent less than 30` = `less than 30`*100 / total, `percent nonelderly` = (100 - `percent elderly`)) %>% 
  left_join(sj_median_income_by_block %>% dplyr::select(-Median_Income)) %>% 
  filter(!is.na(device_count)) 
# get data on vehicles available as vehicles allocation
sj_vehicles_by_block <- getCensus(
    name = "acs/acs5",
    vintage = 2018,
    region = "block group:*", 
    regionin = "state:06+county:085",
    vars = "group(B992512)"
  ) %>% 
  mutate(
    blockgroup =
      paste0(state,county,tract,block_group)
  ) %>% 
  select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>% 
  select(-c(contains("EA"),contains("MA"),contains("M"))) %>% 
  dplyr::select(B992512_001E, blockgroup) %>%
  rename(total_vehicles = B992512_001E, blockgroup = blockgroup) %>%
  left_join(sj_age_by_block %>% dplyr::select_if(!names(.) %in% c("elderly", "percent elderly", "less than 30", "percent less than 30"))) %>%
  mutate(`vehicles per capita` = total_vehicles / total) %>%
  filter(!is.na(device_count)) 

# also get data on vehicles available as households without a vehicle
sj_no_vehicles_by_block <- getCensus(
    name = "acs/acs5",
    vintage = 2018,
    region = "block group:*", 
    regionin = "state:06+county:085",
    vars = "group(B25044)"
  ) %>% 
  mutate(
    blockgroup =
      paste0(state,county,tract,block_group)
  ) %>% 
  select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>% 
  select(-c(contains("EA"),contains("MA"),contains("M"))) %>% 
  gather(key = "variable", value = "estimate", -blockgroup) %>% 
  mutate(label = acs_vars$label[match(variable,acs_vars$name)]) %>% 
  select(-variable) %>%
  separate(label, into = c(NA, NA, NA,"vehicles"), sep = "!!") %>% 
  filter(!is.na(vehicles)) %>%
  group_by(blockgroup, vehicles) %>%
  summarize(grouped_vehicles = sum(estimate)) %>%
  spread(key = vehicles, value = grouped_vehicles) %>%
  mutate(total_nums = `1 vehicle available` + `2 vehicles available` + `3 vehicles available` + `4 vehicles available` + `5 or more vehicles available` + `No vehicle available`, `percent no vehicles` = `No vehicle available`*100 / total_nums, `percent with vehicles` = (100-`percent no vehicles`)) %>%
  left_join(sj_age_by_block %>% dplyr::select_if(!names(.) %in% c("elderly", "percent elderly", "less than 30", "percent less than 30"))) %>%
  filter(!is.na(device_count))
# get data on occupants per room
sj_occupants_per_room_by_block <- getCensus(
    name = "acs/acs5",
    vintage = 2018,
    region = "block group:*", 
    regionin = "state:06+county:085",
    vars = "group(B25014)"
  ) %>% 
  mutate(
    blockgroup =
      paste0(state,county,tract,block_group)
  ) %>% 
  select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>% 
  select(-c(contains("EA"),contains("MA"),contains("M"))) %>% 
  gather(key = "variable", value = "estimate", -blockgroup) %>% 
  mutate(label = acs_vars$label[match(variable,acs_vars$name)]) %>% 
  select(-variable) %>% 
  separate(label, into = c(NA, NA, NA,"occupants per room"), sep = "!!") %>% 
  filter(!is.na(`occupants per room`)) %>%
  group_by(blockgroup, `occupants per room`) %>%
  summarize(estimate_tot = sum(estimate)) %>% 
  spread(key = `occupants per room`, value = estimate_tot) %>%
  mutate(total_nums = `0.50 or less occupants per room` + `0.51 to 1.00 occupants per room` + `1.01 to 1.50 occupants per room` + `1.51 to 2.00 occupants per room` + `2.01 or more occupants per room`, `percent 1 or more` = (`1.01 to 1.50 occupants per room` + `1.51 to 2.00 occupants per room` + `2.01 or more occupants per room`) * 100/ total_nums, `percent less than 1` = (100-`percent 1 or more`)) %>%
  left_join(sj_age_by_block %>% dplyr::select_if(!names(.) %in% c("elderly", "percent elderly", "less than 30", "percent less than 30"))) %>%
  filter(!is.na(device_count)) 

Testing correlations

In the plots below, I show the selected variables against percent of devices completely at home since the shelter-in-place order started, as well as against percent of devices pre-shelter-in-place for comparison.

Age:

# age
sj_age_by_block %>%
  ggplot(aes(
  x = `percent less than 30`,
  y = `% not completely at home`
)) + geom_point() + 
  geom_smooth(method=lm) + 
labs(
    x = "Percent of residents younger than 30",
    y = "Percent devices leaving home on weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Young Age Groups"
  )

young_model <- lm(sj_age_by_block$`% not completely at home` ~ sj_age_by_block$`percent less than 30`)
summary(young_model)
## 
## Call:
## lm(formula = sj_age_by_block$`% not completely at home` ~ sj_age_by_block$`percent less than 30`)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.159  -4.833  -0.233   4.473  39.272 
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                            42.52788    1.52328  27.919  < 2e-16 ***
## sj_age_by_block$`percent less than 30`  0.20725    0.03861   5.367 1.17e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.352 on 567 degrees of freedom
## Multiple R-squared:  0.04835,    Adjusted R-squared:  0.04667 
## F-statistic: 28.81 on 1 and 567 DF,  p-value: 1.167e-07
sj_age_by_block %>% filter(`percent elderly` < 50) %>% # get rid of extreme outliers
  ggplot(aes(
  x = `percent elderly`,
  y = `% not completely at home`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of residents 65 and older",
    y = "Percent devices leaving home on weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Elderly Population"
  )

elderly_model <- lm(`% not completely at home` ~ `percent elderly`, sj_age_by_block %>% filter(`percent elderly` < 50))
summary(elderly_model)
## 
## Call:
## lm(formula = `% not completely at home` ~ `percent elderly`, 
##     data = sj_age_by_block %>% filter(`percent elderly` < 50))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.119  -5.108  -0.259   4.680  36.632 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       52.97809    0.79188   66.90  < 2e-16 ***
## `percent elderly` -0.19522    0.05468   -3.57 0.000387 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.456 on 564 degrees of freedom
## Multiple R-squared:  0.0221, Adjusted R-squared:  0.02037 
## F-statistic: 12.75 on 1 and 564 DF,  p-value: 0.0003869
# compare this to pre-shelter-in-place behavior
sj_age_by_block %>%
  ggplot(aes(
  x = `percent less than 30`,
  y = `% not completely at home pre shelter`
)) + geom_point() + 
  geom_smooth(method=lm) + 
labs(
    x = "Percent of residents younger than 30",
    y = "Percent devices leaving home pre-shelter-in-place",
    title = "San Jose: Staying at Home and Young Age Groups Pre Shelter-in-Place"
  )

young_model2 <- lm(sj_age_by_block$`% not completely at home pre shelter` ~ sj_age_by_block$`percent less than 30`)
summary(young_model2)
## 
## Call:
## lm(formula = sj_age_by_block$`% not completely at home pre shelter` ~ 
##     sj_age_by_block$`percent less than 30`)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -27.0472  -2.4524  -0.1802   2.9456  15.9202 
## 
## Coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                            82.50158    0.77360 106.647  < 2e-16 ***
## sj_age_by_block$`percent less than 30` -0.09199    0.01961  -4.691 3.41e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.242 on 567 degrees of freedom
## Multiple R-squared:  0.03736,    Adjusted R-squared:  0.03566 
## F-statistic: 22.01 on 1 and 567 DF,  p-value: 3.41e-06
sj_age_by_block %>% filter(`percent elderly` < 50) %>% # get rid of extreme outliers
  ggplot(aes(
  x = `percent elderly`,
  y = `% not completely at home pre shelter`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of residents 65 and older",
    y = "Percent devices leaving home on weekdays pre-shelter-in-place",
    title = "San Jose: Staying at Home and Elderly Population Pre Shelter-in-Place"
  )

elderly_model2 <- lm(`% not completely at home pre shelter` ~ `percent elderly`, sj_age_by_block %>% filter(`percent elderly` < 50))
summary(elderly_model2)
## 
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `percent elderly`, 
##     data = sj_age_by_block %>% filter(`percent elderly` < 50))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -26.2336  -2.5173  -0.1708   3.0147  12.2215 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       77.44185    0.39776 194.694  < 2e-16 ***
## `percent elderly`  0.11900    0.02746   4.333 1.74e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.247 on 564 degrees of freedom
## Multiple R-squared:  0.03222,    Adjusted R-squared:  0.0305 
## F-statistic: 18.77 on 1 and 564 DF,  p-value: 1.742e-05

Income:

# income - less than $75000
sj_ami_by_block %>% 
  ggplot(aes(
  x = `% over 75,000`,
  y = `% not completely at home`
)) + geom_point() + 
  geom_smooth(method=lm) +
  labs(
    x = "Percent of housholds with incomes over $75,000 (50% AMI) annually",
    y = "Percent devices leaving home on weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Households Above 50% AMI"
  )

income_75_model <- lm(`% not completely at home` ~ `% over 75,000`, sj_ami_by_block)
summary(income_75_model)
## 
## Call:
## lm(formula = `% not completely at home` ~ `% over 75,000`, data = sj_ami_by_block)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.690  -4.649  -0.541   4.168  34.789 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      64.3419     1.1156   57.68   <2e-16 ***
## `% over 75,000`  -0.2233     0.0172  -12.99   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.447 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.2296, Adjusted R-squared:  0.2282 
## F-statistic: 168.7 on 1 and 566 DF,  p-value: < 2.2e-16
# income - less than $100000
sj_ami_by_block %>% 
  ggplot(aes(
  x = `% over 100,000`,
  y = `% not completely at home`
)) + geom_point() + 
  geom_smooth(method=lm) +
  labs(
    x = "Percent of housholds with incomes over $100,000 (80% AMI) annually",
    y = "Percent devices leaving home on weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Households Below 80% AMI"
  )

income_100_model <- lm(`% not completely at home` ~ `% over 100,000`, sj_ami_by_block)
summary(income_100_model)
## 
## Call:
## lm(formula = `% not completely at home` ~ `% over 100,000`, data = sj_ami_by_block)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.120  -4.507  -0.736   4.007  33.026 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      61.69920    0.86857   71.03   <2e-16 ***
## `% over 100,000` -0.22077    0.01592  -13.87   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.33 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.2536, Adjusted R-squared:  0.2523 
## F-statistic: 192.3 on 1 and 566 DF,  p-value: < 2.2e-16
# income - less than $125000
sj_ami_by_block %>% 
  ggplot(aes(
  x = `% over 125,000`,
  y = `% not completely at home`
)) + geom_point() + 
  geom_smooth(method=lm) +
  labs(
    x = "Percent of housholds with incomes over $125,000 annually",
    y = "Percent devices leaving home on weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Households Below $125,000"
  )

income_125_model <- lm(`% not completely at home` ~ `% over 125,000`, sj_ami_by_block)
summary(income_125_model)
## 
## Call:
## lm(formula = `% not completely at home` ~ `% over 125,000`, data = sj_ami_by_block)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.135  -4.477  -0.883   4.279  31.958 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      59.89204    0.72846   82.22   <2e-16 ***
## `% over 125,000` -0.22989    0.01608  -14.30   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.272 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.2653, Adjusted R-squared:  0.264 
## F-statistic: 204.4 on 1 and 566 DF,  p-value: < 2.2e-16
# compare to pre shelter in place
sj_ami_by_block %>% 
  ggplot(aes(
  x = `% over 75,000`,
  y = `% not completely at home pre shelter`
)) + geom_point() + 
  geom_smooth(method=lm) +
  labs(
    x = "Percent of housholds with incomes over $75,000 (50% AMI) annually",
    y = "Percent devices leaving home on weekdays pre-shelter-in-place",
    title = "San Jose: Staying at Home and Households Above 50% AMI Pre Shelter-in-Place"
  )

income_75_model2 <- lm(`% not completely at home pre shelter` ~ `% over 75,000`, sj_ami_by_block)
summary(income_75_model2)
## 
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `% over 75,000`, 
##     data = sj_ami_by_block)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -27.3258  -2.4328   0.0752   2.8199  14.2150 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     74.335377   0.615515 120.769  < 2e-16 ***
## `% over 75,000`  0.074344   0.009487   7.836 2.32e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.109 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.09787,    Adjusted R-squared:  0.09627 
## F-statistic:  61.4 on 1 and 566 DF,  p-value: 2.315e-14
# income - less than $100000
sj_ami_by_block %>% 
  ggplot(aes(
  x = `% over 100,000`,
  y = `% not completely at home pre shelter`
)) + geom_point() + 
  geom_smooth(method=lm) +
  labs(
    x = "Percent of housholds with incomes over $100,000 (80% AMI) annually",
    y = "Percent devices leaving home on weekdays pre-shelter-in-place",
    title = "San Jose: Staying Home and Households Below 80% AMI Pre Shelter-in-Place"
  )

income_100_model2 <- lm(`% not completely at home pre shelter` ~ `% over 100,000`, sj_ami_by_block)
summary(income_100_model2)
## 
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `% over 100,000`, 
##     data = sj_ami_by_block)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -27.3698  -2.4317   0.0472   2.7043  14.3737 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      75.026275   0.481081 155.953   <2e-16 ***
## `% over 100,000`  0.077203   0.008818   8.755   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.06 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.1193, Adjusted R-squared:  0.1177 
## F-statistic: 76.66 on 1 and 566 DF,  p-value: < 2.2e-16
# over 125000
sj_ami_by_block %>% 
  ggplot(aes(
  x = `% over 125,000`,
  y = `% not completely at home pre shelter`
)) + geom_point() + 
  geom_smooth(method=lm) +
  labs(
    x = "Percent of housholds with incomes over $125,000 annually",
    y = "Percent devices leaving home on weekdays pre-shelter-in-place",
    title = "San Jose: Social Distancing and Households Below $125,000 Pre Shelter-in-Place"
  )

income_125_model2 <- lm(`% not completely at home pre shelter` ~ `% over 125,000`, sj_ami_by_block)
summary(income_125_model2)
## 
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `% over 125,000`, 
##     data = sj_ami_by_block)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -27.2220  -2.3676   0.1555   2.5383  13.9673 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      75.432749   0.401316 187.964   <2e-16 ***
## `% over 125,000`  0.085872   0.008858   9.694   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.006 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.1424, Adjusted R-squared:  0.1409 
## F-statistic: 93.98 on 1 and 566 DF,  p-value: < 2.2e-16

Language:

# language
sj_lang_by_block %>% 
  ggplot(aes(
  x = `% speaking english > well`,
  y = `% not completely at home`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of individuals speaking English well",
    y = "Percent devices leaving home on weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and English Language Ability"
  )

english_ability_model <- lm(`% not completely at home` ~ `% speaking english > well`, sj_lang_by_block)
summary(english_ability_model)
## 
## Call:
## lm(formula = `% not completely at home` ~ `% speaking english > well`, 
##     data = sj_lang_by_block)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.526  -4.999  -0.397   3.914  38.388 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 70.36942    3.36653  20.903  < 2e-16 ***
## `% speaking english > well` -0.22415    0.03775  -5.938 5.02e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.307 on 567 degrees of freedom
## Multiple R-squared:  0.05855,    Adjusted R-squared:  0.05689 
## F-statistic: 35.26 on 1 and 567 DF,  p-value: 5.02e-09
sj_lang_by_block %>% 
  ggplot(aes(
  x = `% not speaking spanish`,
  y = `% not completely at home`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of individuals not speaking Spanish",
    y = "Percent devices leaving home on weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Spanish Language Ability"
  )

spanish_speaking_model <- lm(`% not completely at home` ~ `% not speaking spanish`, sj_lang_by_block)
summary(spanish_speaking_model)
## 
## Call:
## lm(formula = `% not completely at home` ~ `% not speaking spanish`, 
##     data = sj_lang_by_block)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.705  -4.658  -0.741   4.051  37.548 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              63.80076    1.32061   48.31   <2e-16 ***
## `% not speaking spanish` -0.17128    0.01645  -10.41   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.845 on 567 degrees of freedom
## Multiple R-squared:  0.1605, Adjusted R-squared:  0.159 
## F-statistic: 108.4 on 1 and 567 DF,  p-value: < 2.2e-16
# compare to pre shelter in place
sj_lang_by_block %>% 
  ggplot(aes(
  x = `% speaking english > well`,
  y = `% not completely at home pre shelter`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of individuals speaking English well",
    y = "Percent devices leaving home on weekdays pre-shelter-in-place",
    title = "San Jose: Staying at Home and English Language Ability Pre Shelter-in-Place"
  )

english_ability_model2 <- lm(`% not completely at home pre shelter` ~ `% speaking english > well`, sj_lang_by_block)
summary(english_ability_model2)
## 
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `% speaking english > well`, 
##     data = sj_lang_by_block)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -27.7048  -2.2760   0.0598   2.8560  10.3112 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 64.51050    1.64215  39.284   <2e-16 ***
## `% speaking english > well`  0.16300    0.01841   8.853   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.052 on 567 degrees of freedom
## Multiple R-squared:  0.1214, Adjusted R-squared:  0.1199 
## F-statistic: 78.37 on 1 and 567 DF,  p-value: < 2.2e-16
sj_lang_by_block %>% 
  ggplot(aes(
  x = `% not speaking spanish`,
  y = `% not completely at home pre shelter`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of individuals not speaking Spanish",
    y = "Percent devices leaving home on weekdays pre shelter-in-place",
    title = "San Jose: Staying at Home and Spanish Language Ability Pre Shelter-in-Place"
  )

spanish_speaking_model2 <- lm(`% not completely at home pre shelter` ~ `% not speaking spanish`, sj_lang_by_block)
summary(spanish_speaking_model2)
## 
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `% not speaking spanish`, 
##     data = sj_lang_by_block)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -27.6257  -2.4274   0.0822   2.7537  10.4706 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              73.221669   0.683770  107.08   <2e-16 ***
## `% not speaking spanish`  0.073935   0.008518    8.68   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.062 on 567 degrees of freedom
## Multiple R-squared:  0.1173, Adjusted R-squared:  0.1157 
## F-statistic: 75.34 on 1 and 567 DF,  p-value: < 2.2e-16

Occupants per room:

# occupants per room
sj_occupants_per_room_by_block %>% 
  ggplot(aes(
  x = `percent less than 1`,
  y = `% not completely at home`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of households with 1 or fewer occupant per room",
    y = "Percent devices leaving home on weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Room Occupancy"
  )

occupants_model <- lm(`% not completely at home` ~ `percent less than 1`, sj_occupants_per_room_by_block)
summary(occupants_model)
## 
## Call:
## lm(formula = `% not completely at home` ~ `percent less than 1`, 
##     data = sj_occupants_per_room_by_block)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.842  -4.849  -0.286   4.137  34.977 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           71.54113    2.97589   24.04  < 2e-16 ***
## `percent less than 1` -0.23399    0.03277   -7.14 2.89e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.126 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.08262,    Adjusted R-squared:  0.081 
## F-statistic: 50.97 on 1 and 566 DF,  p-value: 2.885e-12
# compare to pre shelter in place
sj_occupants_per_room_by_block %>% 
  ggplot(aes(
  x = `percent less than 1`,
  y = `% not completely at home pre shelter`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of households with 1 or fewer occupant per room",
    y = "Percent devices leaving home on weekdays pre shelter-in-place",
    title = "San Jose: Staying at Home and Room Occupancy Pre Shelter-in-Place"
  )

occupants_model2 <- lm(`% not completely at home pre shelter` ~ `percent less than 1`, sj_occupants_per_room_by_block)
summary(occupants_model2)
## 
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `percent less than 1`, 
##     data = sj_occupants_per_room_by_block)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -27.1950  -2.5198  -0.0566   2.8717  14.3938 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           65.90546    1.48474  44.389   <2e-16 ***
## `percent less than 1`  0.14478    0.01635   8.854   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.054 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.1217, Adjusted R-squared:  0.1201 
## F-statistic:  78.4 on 1 and 566 DF,  p-value: < 2.2e-16

Vehicle ownership:

# vehicles
sj_vehicles_by_block %>% 
  ggplot(aes(
  x = `vehicles per capita`,
  y = `% not completely at home`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Number of vehicles per capita",
    y = "Percent devices leaving home on weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Vehicles Per Capita"
  )

# vehicles - percent with no vehicles
sj_no_vehicles_by_block %>% 
  ggplot(aes(
  x = `percent with vehicles`,
  y = `% not completely at home`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of housholds with vehicles available",
    y = "Percent devices leaving home on weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Vehicle Availability"
  )

vehicles_model <- lm(`% not completely at home` ~ `percent with vehicles`, sj_no_vehicles_by_block)
summary(vehicles_model)
## 
## Call:
## lm(formula = `% not completely at home` ~ `percent with vehicles`, 
##     data = sj_no_vehicles_by_block)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.696  -5.096  -0.235   4.625  38.404 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             78.75974    5.18753  15.183  < 2e-16 ***
## `percent with vehicles` -0.29763    0.05439  -5.473 6.67e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.268 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.05025,    Adjusted R-squared:  0.04858 
## F-statistic: 29.95 on 1 and 566 DF,  p-value: 6.672e-08
# compare to pre shelter in place
sj_no_vehicles_by_block %>% 
  ggplot(aes(
  x = `percent with vehicles`,
  y = `% not completely at home pre shelter`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of housholds with vehicles available",
    y = "Percent devices leaving home on weekdays pre shelter-in-place",
    title = "San Jose: Social Distancing and Vehicle Availability Pre Shelter-in-Place"
  )

vehicles_model2 <- lm(`% not completely at home pre shelter` ~ `percent with vehicles`, sj_no_vehicles_by_block)
summary(vehicles_model2)
## 
## Call:
## lm(formula = `% not completely at home pre shelter` ~ `percent with vehicles`, 
##     data = sj_no_vehicles_by_block)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -24.6695  -2.7766   0.0556   2.9241  10.8094 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             65.70738    2.65606  24.739  < 2e-16 ***
## `percent with vehicles`  0.13931    0.02785   5.003 7.55e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.233 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.04235,    Adjusted R-squared:  0.04066 
## F-statistic: 25.03 on 1 and 566 DF,  p-value: 7.55e-07

Multiple regression analysis with income, age, language, and occupants per room

# multiple regression 
modeltest <- lm(sj_ami_by_block$`% not completely at home` ~ sj_ami_by_block$`% over 125,000` + sj_age_by_block$`percent less than 30` + sj_lang_by_block$`% speaking english > well` + sj_occupants_per_room_by_block$`percent less than 1`)
summary(modeltest)
## 
## Call:
## lm(formula = sj_ami_by_block$`% not completely at home` ~ sj_ami_by_block$`% over 125,000` + 
##     sj_age_by_block$`percent less than 30` + sj_lang_by_block$`% speaking english > well` + 
##     sj_occupants_per_room_by_block$`percent less than 1`)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.884  -4.568  -0.779   4.270  31.930 
## 
## Coefficients:
##                                                       Estimate Std. Error
## (Intercept)                                          50.772738   4.743391
## sj_ami_by_block$`% over 125,000`                     -0.250139   0.021391
## sj_age_by_block$`percent less than 30`                0.041453   0.041836
## sj_lang_by_block$`% speaking english > well`          0.090269   0.045832
## sj_occupants_per_room_by_block$`percent less than 1`  0.003916   0.045858
##                                                      t value Pr(>|t|)    
## (Intercept)                                           10.704   <2e-16 ***
## sj_ami_by_block$`% over 125,000`                     -11.694   <2e-16 ***
## sj_age_by_block$`percent less than 30`                 0.991   0.3222    
## sj_lang_by_block$`% speaking english > well`           1.970   0.0494 *  
## sj_occupants_per_room_by_block$`percent less than 1`   0.085   0.9320    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.253 on 563 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.273,  Adjusted R-squared:  0.2679 
## F-statistic: 52.86 on 4 and 563 DF,  p-value: < 2.2e-16

Testing animating the plot

# collect the demographic variables
sj_dem_distancing <- sj_internet_by_block %>% 
  dplyr::select(`percent high speed`, `% not completely at home`, blockgroup) %>% 
  left_join(sj_education_by_block %>% dplyr::select(blockgroup, `percent associates or higher`)) %>% 
  left_join(sj_ami_by_block %>% dplyr::select(blockgroup, `% over 125,000`)) %>% 
  left_join(sj_age_by_block %>% dplyr::select(blockgroup, `percent less than 30`)) %>% 
  left_join(sj_lang_by_block %>% dplyr::select(blockgroup, `% not speaking spanish`)) %>% 
  left_join(sj_no_vehicles_by_block %>% dplyr::select(blockgroup, `percent with vehicles`)) %>%
  left_join(sj_occupants_per_room_by_block %>% dplyr::select(blockgroup, `percent less than 1`))

# another collection for pre shelter in place behavior
sj_dem_distancing_pre_shelter <- sj_dem_distancing %>% 
  dplyr::select(-`% not completely at home`) %>%
  left_join(sj_internet_by_block %>% dplyr::select(`% not completely at home pre shelter`, blockgroup))

# relabel column for leaving home
colnames(sj_dem_distancing_pre_shelter)[9] <- "% not completely at home"

sj_dem_distancing[is.na(sj_dem_distancing)] <- 0
sj_dem_distancing_pre_shelter[is.na(sj_dem_distancing_pre_shelter)] <- 0

# add column indicating before or after shelter in place, then bind the two sets of data
sj_dem_distancing_pre_shelter <- sj_dem_distancing_pre_shelter %>% 
  mutate(
    income_trendline =
      fitted(lm(sj_dem_distancing_pre_shelter$`% not completely at home` ~ sj_dem_distancing_pre_shelter$`% over 125,000`))
  ) %>% 
  cbind(`Before or After Shelter-in-Place` = "before")
sj_dem_distancing <-
  sj_dem_distancing %>%
  mutate(
    income_trendline =
      fitted(lm(sj_dem_distancing$`% not completely at home` ~ sj_dem_distancing$`% over 125,000`))
  ) %>% 
  cbind(`Before or After Shelter-in-Place` = "after") %>% 
  rbind(sj_dem_distancing_pre_shelter)

# try animating
fig <- 
  plot_ly (sj_dem_distancing) %>%
    add_trace(
      x = ~`% over 125,000`, 
      y = ~`% not completely at home`, 
      frame = ~`Before or After Shelter-in-Place`, 
      type = 'scatter', 
      mode = 'markers'
    ) %>% 
    add_trace(
      x = ~`% over 125,000`,
      y = ~income_trendline,
      type = 'scatter',
      mode = 'lines',
      line = list(size = 5, color = 'rgba(255, 165, 0, 1)'),
      frame = ~`Before or After Shelter-in-Place`
    ) %>% 
  animation_button(visible = F)
fig
# # save as rds
# saveRDS(sj_dem_distancing, "/Users/simonespeizer/pCloud Drive/Shared/SFBI/Restricted Data Library/Safegraph/covid19analysis/sj_sd_dem_data.rds")


# fig <- plot_ly(sj_dem_distancing) %>% 
#   add_trace(
#     x = ~`% over 125,000`,
#     y = ~`% not completely at home`,
#     frame = ~`Before or After Shelter-in-Place`,
#     type = "scatter",
#     mode = "markers",
#     name = "Under $125,000",
#     marker = list(size = 5, color = 'rgba(50, 150, 200, 1)'),
#     visible = T
#   ) %>% 
#   add_trace(
#     x = ~`% over 125,000`,
#     y = fitted(lm(sj_dem_distancing$`% not completely at home` ~ sj_dem_distancing$`% over 125,000`)),
#     name = 'trendline',
#     mode = 'lines',
#     line = list(size = 5, color = 'rgba(255, 165, 0, 1)'),
#     frame = ~`Before or After Shelter-in-Place`,
#     visible = F
#   ) %>%
#   add_trace(
#     x = ~`% not speaking spanish`,
#     y = ~`% not completely at home`,
#     frame = ~`Before or After Shelter-in-Place`,
#     name = "speak spanish",
#     marker = list(size = 5, color = 'rgba(50, 150, 200, 1)'),
#     visible = F
#   ) %>% 
#   add_trace(
#     x = ~`% not speaking spanish`,
#     y = fitted(lm(sj_dem_distancing$`% not completely at home` ~ sj_dem_distancing$`% not speaking spanish`)),
#     name = 'trendline',
#     mode = 'lines',
#     line = list(size = 5, color = 'rgba(255, 165, 0, 1)'),
#     frame = ~`Before or After Shelter-in-Place`,
#     visible = F
#   ) %>% 
#   add_trace(
#     x = ~`percent associates or higher`,
#     y = ~`% not completely at home`,
#     frame = ~`Before or After Shelter-in-Place`,
#     name = "percent higher degree",
#     marker = list(size = 5, color = 'rgba(50, 150, 200, 1)'),
#     visible = F
#   ) %>% 
#   add_trace(
#     x = ~`percent associates or higher`,
#     y = fitted(lm(sj_dem_distancing$`% not completely at home` ~ sj_dem_distancing$`percent associates or higher`)),
#     name = 'trendline',
#     mode = 'lines',
#     line = list(size = 5, color = 'rgba(255, 165, 0, 1)'),
#     frame = ~`Before or After Shelter-in-Place`,
#     visible = F
#   ) %>%
#   add_trace(
#     x = ~`percent high speed`,
#     y = ~`% not completely at home`,
#     frame = ~`Before or After Shelter-in-Place`,
#     name = "percent high speed internet access",
#     marker = list(size = 5, color = 'rgba(50, 150, 200, 1)'),
#     visible = F
#   ) %>% 
#   add_trace(
#     x = ~`percent high speed`,
#     y = fitted(lm(sj_dem_distancing$`% not completely at home` ~ sj_dem_distancing$`percent high speed`)),
#     name = 'trendline',
#     mode = 'lines',
#     line = list(size = 5, color = 'rgba(255, 165, 0, 1)'),
#     frame = ~sj_dem_distancing$`Before or After Shelter-in-Place`,
#     visible = F
#   ) %>%
#   add_trace(
#     x = ~`percent less than 30`,
#     y = ~`% not completely at home`,
#     frame = ~`Before or After Shelter-in-Place`,
#     name = "percent less than 30",
#     marker = list(size = 5, color = 'rgba(50, 150, 200, 1)'),
#     visible = F
#   ) %>% 
#   add_trace(
#     x = ~`percent less than 30`,
#     y = fitted(lm(sj_dem_distancing$`% not completely at home` ~ sj_dem_distancing$`percent less than 30`)),
#     name = 'trendline',
#     mode = 'lines',
#     line = list(size = 5, color = 'rgba(255, 165, 0, 1)'),
#     frame = ~`Before or After Shelter-in-Place`,
#     visible = F
#   ) %>%
#   layout(
#     updatemenus = list(
#       list(
#         active = 0,
#         type = 'buttons',
#         buttons = list(
#           list(
#             label = "Households Under $125,000",
#             method = 'update',
#             args = list(list(visible = c(T, T, F, F, F, F, F, F, F, F)),
#                         list(title = "Under $125,000",
#                              xaxis = list(title = "% Households Under $125,000 in Income")))),
#           list(
#             label = "Speaking Spanish",
#             method = 'update',
#             args = list(list(visible = c(F, F, T, T, F, F, F, F, F, F)),
#                         list(title = "Not Speaking Spanish",
#                              xaxis = list(title = "% Residents Not Speaking Spanish")))),
#           list(
#             label = "Education Level",
#             method = 'update',
#             args= list(list(visible = c(F, F, F, F, T, T, F, F, F, F)),
#                        list(xaxis = list(title = "% Residents With Associate's Degree or Higher")))),
#           list(
#             label = "High Speed Internet",
#             method = 'update',
#             args= list(list(visible = c(F, F, F, F, F, F, T, T, F, F)),
#                        list(xaxis = list(title = "% Households With High Speed Internet Access")))),
#           list(
#             label = "Young Population",
#             method = 'update',
#             args= list(list(visible = c(F, F, F, F, F, F, T, T, F, F)),
#                        list(xaxis = list(title = "% Residents Under Age 30"))))
#           )
#           )
#         ),
#     yaxis = list(title = "% Residents Leaving Home", 
#                  font = list(size = 15)),
#     showlegend = FALSE
#       )
# fig

Experimentation

Experimentation with other variables and other ways of analyzing the social distancing data. First I look at a few other possible variables. I start with units in the structure.

# try getting other variables
# get data on units in structure
sj_units_in_structure_by_block <- getCensus(
    name = "acs/acs5",
    vintage = 2018,
    region = "block group:*", 
    regionin = "state:06+county:085",
    vars = "group(B25024)"
  ) %>% 
  mutate(
    blockgroup =
      paste0(state,county,tract,block_group)
  ) %>% 
  select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>% 
  select(-c(contains("EA"),contains("MA"),contains("M"))) %>% 
  gather(key = "variable", value = "estimate", -blockgroup) %>% 
  mutate(label = acs_vars$label[match(variable,acs_vars$name)]) %>% 
  select(-variable) %>% 
  separate(label, into = c(NA, NA, "units"), sep = "!!") %>% 
  filter(!is.na(units)) %>%
  spread(key = units, value = estimate) %>%
  mutate(total_nums = `1, attached` + `1, detached` + `10 to 19` + `2` + `20 to 49`+ `3 or 4` + `5 to 9`+ `50 or more`+ `Boat, RV, van, etc.`+ `Mobile home`, `percent 20 or more` = (`20 to 49`+`50 or more`)* 100/ total_nums, `percent 1 only` = (`1, attached` + `1, detached`)*100/total_nums, `percent > 1` = 100 - `percent 1 only`) %>%
  left_join(sj_age_by_block %>% dplyr::select_if(!names(.) %in% c("elderly", "percent elderly", "less than 30", "percent less than 30"))) %>%
  filter(!is.na(device_count))

# plot 
sj_units_in_structure_by_block %>% 
  ggplot(aes(
  x = `percent 20 or more`,
  y = `% not completely at home`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of structures with more than 20 units",
    y = "Percent devices leaving home on weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and 20 or More Units Per Structure"
  )

summary(lm(`% not completely at home` ~ `percent 20 or more`, sj_units_in_structure_by_block))
## 
## Call:
## lm(formula = `% not completely at home` ~ `percent 20 or more`, 
##     data = sj_units_in_structure_by_block)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.765  -5.190  -0.201   4.843  37.335 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          50.06531    0.40940 122.290   <2e-16 ***
## `percent 20 or more`  0.03712    0.02052   1.809    0.071 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.459 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.005748,   Adjusted R-squared:  0.003992 
## F-statistic: 3.272 on 1 and 566 DF,  p-value: 0.07099
sj_units_in_structure_by_block %>% 
  ggplot(aes(
  x = `percent 1 only`,
  y = `% not completely at home`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of structures with only one unit",
    y = "Percent devices leaving home on weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Only 1 Unit Per Structure"
  )

summary(lm(`% not completely at home` ~ `percent 1 only`, sj_units_in_structure_by_block))
## 
## Call:
## lm(formula = `% not completely at home` ~ `percent 1 only`, data = sj_units_in_structure_by_block)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.608  -5.078  -0.224   4.405  38.261 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      55.24549    0.88920  62.130  < 2e-16 ***
## `percent 1 only` -0.06648    0.01132  -5.872 7.33e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.237 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.05743,    Adjusted R-squared:  0.05576 
## F-statistic: 34.48 on 1 and 566 DF,  p-value: 7.328e-09

Household type and size:

# load data on household type and size
sj_house_size_type_by_block <- getCensus(
    name = "acs/acs5",
    vintage = 2018,
    region = "block group:*", 
    regionin = "state:06+county:085",
    vars = "group(B11016)"
  ) %>% 
  mutate(
    blockgroup =
      paste0(state,county,tract,block_group)
  ) %>% 
  select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>% 
  select(-c(contains("EA"),contains("MA"),contains("M"))) %>% 
  gather(key = "variable", value = "estimate", -blockgroup) %>% 
  mutate(label = acs_vars$label[match(variable,acs_vars$name)]) %>% 
  select(-variable) %>% 
  separate(label, into = c(NA, NA, "type", "size"), sep = "!!") %>% 
  filter(!is.na(type))


# household type
sj_house_type_by_block <- sj_house_size_type_by_block %>% 
  filter(is.na(size)) %>% 
  dplyr::select(-size) %>%
  spread(key = type, value = estimate) %>% 
  mutate(`total households` = `Family households` + `Nonfamily households`, `percent nonfamily` = `Nonfamily households` / `total households`) %>%
  left_join(sj_age_by_block %>% dplyr::select_if(!names(.) %in% c("elderly", "percent elderly", "less than 30", "percent less than 30"))) %>%
  filter(!is.na(device_count))

sj_house_type_by_block %>% 
  ggplot(aes(
  x = `percent nonfamily`,
  y = `% not completely at home`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent nonfamily households",
    y = "Percent devices leaving home on weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Household Type"
  )

summary(lm(`% not completely at home` ~ `percent nonfamily`, sj_house_type_by_block))
## 
## Call:
## lm(formula = `% not completely at home` ~ `percent nonfamily`, 
##     data = sj_house_type_by_block)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.828  -5.089   0.052   4.597  38.208 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          47.7692     0.6403  74.609  < 2e-16 ***
## `percent nonfamily`  11.0270     2.2223   4.962 9.24e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.305 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.04169,    Adjusted R-squared:  0.03999 
## F-statistic: 24.62 on 1 and 566 DF,  p-value: 9.245e-07
# household size
sj_house_size_by_block <- sj_house_size_type_by_block %>% 
  filter(!is.na(size)) %>% 
  dplyr::select(-type) %>%
  group_by(blockgroup, size) %>%
  summarize(`total of this size` = sum(estimate)) %>% 
  spread(key = size, value = `total of this size`) %>%
  mutate(total_nums = `1-person household` + `2-person household` + `3-person household` + `4-person household` + `5-person household`+ `6-person household` + `7-or-more person household`, `percent 5 or more` = (`5-person household`+`6-person household` + `7-or-more person household`)* 100/ total_nums, `percent 1 or 2 only` = (`1-person household` + `2-person household`)*100/total_nums) %>%
  left_join(sj_age_by_block %>% dplyr::select_if(!names(.) %in% c("elderly", "percent elderly", "less than 30", "percent less than 30"))) %>%
  filter(!is.na(device_count))

sj_house_size_by_block %>% 
  ggplot(aes(
  x = `percent 5 or more`,
  y = `% not completely at home`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of households with 5 or more people",
    y = "Percent devices leaving home on weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Households With 5 or More"
  )

summary(lm(`% not completely at home` ~ `percent 5 or more`, sj_house_size_by_block))
## 
## Call:
## lm(formula = `% not completely at home` ~ `percent 5 or more`, 
##     data = sj_house_size_by_block)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.419  -4.958  -0.605   4.419  37.588 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         48.67686    0.56617  85.976  < 2e-16 ***
## `percent 5 or more`  0.10054    0.02541   3.957 8.55e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.369 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.02692,    Adjusted R-squared:  0.0252 
## F-statistic: 15.66 on 1 and 566 DF,  p-value: 8.545e-05
sj_house_size_by_block %>% 
  ggplot(aes(
  x = `percent 1 or 2 only`,
  y = `% not completely at home`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of households with 1 or 2 people",
    y = "Percent devices leaving home on weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Small Household Size"
  )

summary(lm(`% not completely at home` ~ `percent 1 or 2 only`, sj_house_size_by_block))
## 
## Call:
## lm(formula = `% not completely at home` ~ `percent 1 or 2 only`, 
##     data = sj_house_size_by_block)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.408  -5.383  -0.223   4.811  37.673 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           49.45494    0.98251  50.335   <2e-16 ***
## `percent 1 or 2 only`  0.02185    0.02043   1.069    0.285    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.475 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.002016,   Adjusted R-squared:  0.0002529 
## F-statistic: 1.143 on 1 and 566 DF,  p-value: 0.2854

Next I consider different ways of looking at the social distancing data. First I try distance traveled.

# try other ways of looking at the social distancing data
# first look at total distance traveled
sj_sd_distance <- sj_socialdistancing %>% 
  filter(date > shelter_start) %>% 
  group_by(origin_census_block_group) %>% 
  summarize(total_dist_traveled = sum(distance_traveled_from_home), device_count = sum(device_count)) %>%
  mutate(total_dist_per_device = total_dist_traveled / device_count)

sj_distance_testing <- left_join(sj_ami_by_block, sj_sd_distance, by = c("blockgroup" = "origin_census_block_group")) %>% left_join(sj_age_by_block %>% select(blockgroup, `percent less than 30`))

sj_distance_testing %>% filter(total_dist_per_device < 500)  %>% 
  ggplot(aes(
  x = `% over 75,000`,
  y = total_dist_per_device
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of housholds with incomes over $75,000 (50% AMI) annually",
    y = "Average distance traveled per device during weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Income, Distance Metric"
  )

This is very skewed by outliers, and probably not a useful metric.

Now I consider including devices that traveled <1km as staying at (or near) home.

sj_sd_range <- sj_socialdistancing %>% 
  filter(weekend == F) %>% 
  filter(date > shelter_start) %>%
  mutate(travel_buckets_split = lapply(bucketed_distance_traveled, function(x) strsplit(x, "<1000")[[1]][2]), less_than_1km = lapply(travel_buckets_split, function(x) strsplit(x, ":")[[1]][2]), less_than_1km = lapply(less_than_1km, function(x) strsplit(x, ",")[[1]][1])) %>%
  mutate(less_than_1km = lapply(less_than_1km, function(x) str_remove(x, "[}]")))  %>% # clean a bit more
  mutate(less_than_1km = as.numeric(less_than_1km), less_than_1km = replace_na(less_than_1km, 0)) %>% 
  mutate(home_or_1km = completely_home_device_count + less_than_1km) %>% 
  group_by(origin_census_block_group) %>% 
  summarize(home_or_1km = sum(home_or_1km), device_count = sum(device_count)) %>% 
  mutate(`% Within 1km of Home` = (home_or_1km/device_count*100) %>% round(1), `% farther than 1km` = (100-`% Within 1km of Home`))

# join this with other data
sj_1km_testing <- left_join(sj_ami_by_block, sj_sd_range, by = c("blockgroup" = "origin_census_block_group")) %>% 
  left_join(sj_occupants_per_room_by_block %>% dplyr::select(`percent less than 1`, blockgroup)) %>%
  left_join(sj_age_by_block %>% dplyr::select(`percent less than 30`, blockgroup)) %>%
  left_join(sj_lang_by_block %>% dplyr::select(`% speaking english > well`, blockgroup)) 

# plot with income
sj_1km_testing %>%  
  ggplot(aes(
  x = `% over 75,000`,
  y = `% farther than 1km`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of housholds with incomes over $75,000 (50% AMI) annually",
    y = "% of devices going farther than 1km of home, weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Income, 1km Range"
  )

summary(lm(`% farther than 1km` ~ `% over 75,000`, sj_1km_testing))
## 
## Call:
## lm(formula = `% farther than 1km` ~ `% over 75,000`, data = sj_1km_testing)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.159  -4.993  -0.401   4.266  41.458 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      56.2525     1.1743   47.90   <2e-16 ***
## `% over 75,000`  -0.2343     0.0181  -12.95   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.839 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.2285, Adjusted R-squared:  0.2271 
## F-statistic: 167.6 on 1 and 566 DF,  p-value: < 2.2e-16
# plot with age
sj_1km_testing %>%  
  ggplot(aes(
  x = `percent less than 30`,
  y = `% farther than 1km`
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of people younger than 30",
    y = "Percent of devices farther than 1km of home during weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Age, 1km Range"
  )

summary(lm(`% farther than 1km` ~ `percent less than 30`, sj_1km_testing))
## 
## Call:
## lm(formula = `% farther than 1km` ~ `percent less than 30`, data = sj_1km_testing)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -25.413  -5.309  -0.304   4.766  42.271 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            31.21800    1.57767  19.787  < 2e-16 ***
## `percent less than 30`  0.27326    0.03999   6.833 2.15e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.651 on 567 degrees of freedom
## Multiple R-squared:  0.07608,    Adjusted R-squared:  0.07445 
## F-statistic: 46.69 on 1 and 567 DF,  p-value: 2.154e-11
# run multiple regression model
modeltest2 <- lm(sj_1km_testing$`% farther than 1km` ~ sj_1km_testing$`% over 75,000` + sj_1km_testing$`percent less than 30` + sj_1km_testing$`% speaking english > well` + sj_1km_testing$`percent less than 1`)
summary(modeltest2)
## 
## Call:
## lm(formula = sj_1km_testing$`% farther than 1km` ~ sj_1km_testing$`% over 75,000` + 
##     sj_1km_testing$`percent less than 30` + sj_1km_testing$`% speaking english > well` + 
##     sj_1km_testing$`percent less than 1`)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.400  -4.765  -0.793   4.599  40.972 
## 
## Coefficients:
##                                            Estimate Std. Error t value Pr(>|t|)
## (Intercept)                                50.12037    4.91129  10.205  < 2e-16
## sj_1km_testing$`% over 75,000`             -0.20967    0.02267  -9.249  < 2e-16
## sj_1km_testing$`percent less than 30`       0.14057    0.04466   3.148  0.00173
## sj_1km_testing$`% speaking english > well` -0.02227    0.04822  -0.462  0.64432
## sj_1km_testing$`percent less than 1`        0.01303    0.04897   0.266  0.79024
##                                               
## (Intercept)                                ***
## sj_1km_testing$`% over 75,000`             ***
## sj_1km_testing$`percent less than 30`      ** 
## sj_1km_testing$`% speaking english > well`    
## sj_1km_testing$`percent less than 1`          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.77 on 563 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.246,  Adjusted R-squared:  0.2406 
## F-statistic: 45.92 on 4 and 563 DF,  p-value: < 2.2e-16

It looks like the fit of these selected variables is slightly better for the social distancing data based on not traveling farther than 1km.

Now I also consider “non-work” behavior.

sj_nonworking_by_block <- sj_socialdistancing %>% 
  filter(weekend == F) %>% 
  filter(date > shelter_start) %>%
  mutate(nonworking = device_count - completely_home_device_count - part_time_work_behavior_devices - full_time_work_behavior_devices) %>%
  group_by(origin_census_block_group) %>%
  summarize(nonworking_count = sum(nonworking), total_device = sum(device_count)) %>% 
  mutate(nonworking_percent = nonworking_count*100 / total_device, percent_only_work_home = 100-nonworking_percent) %>%
  left_join(sj_1km_testing %>% dplyr::select(`% over 75,000`, `percent less than 30`, `% speaking english > well`, `percent less than 1`, blockgroup), by = c("origin_census_block_group" = "blockgroup"))


# plot against age and income
sj_nonworking_by_block %>%  
  ggplot(aes(
  x = `% over 75,000`,
  y = nonworking_percent
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of housholds with incomes over $75,000 (50% AMI) annually",
    y = "Percent of devices leaving home for non-work purposes during weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Income, Nonworking Behavior"
  )

summary(lm(nonworking_percent ~ `% over 75,000`, sj_nonworking_by_block))
## 
## Call:
## lm(formula = nonworking_percent ~ `% over 75,000`, data = sj_nonworking_by_block)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -22.1765  -3.1343  -0.1581   3.1028  17.2184 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     41.24081    0.77032   53.54   <2e-16 ***
## `% over 75,000` -0.12775    0.01187  -10.76   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.142 on 566 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.1698, Adjusted R-squared:  0.1683 
## F-statistic: 115.8 on 1 and 566 DF,  p-value: < 2.2e-16
sj_nonworking_by_block %>%  
  ggplot(aes(
  x = `percent less than 30`,
  y = nonworking_percent
)) + geom_point() + 
  geom_smooth(method=lm) + 
  labs(
    x = "Percent of people younger than 30",
    y = "Percent of devices leaving home for non-work purposes during weekdays since shelter-in-place",
    title = "San Jose: Social Distancing and Age, Nonworking Behavior"
  )

summary(lm(nonworking_percent ~ `percent less than 30`, sj_nonworking_by_block))
## 
## Call:
## lm(formula = nonworking_percent ~ `percent less than 30`, data = sj_nonworking_by_block)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -22.6311  -3.5775  -0.2204   3.4844  17.7223 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            31.16288    1.02739  30.332   <2e-16 ***
## `percent less than 30`  0.05479    0.02604   2.104   0.0358 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.633 on 567 degrees of freedom
## Multiple R-squared:  0.007747,   Adjusted R-squared:  0.005997 
## F-statistic: 4.427 on 1 and 567 DF,  p-value: 0.03582
# multiple regression model
modeltest3 <- lm(sj_nonworking_by_block$nonworking_percent ~ sj_nonworking_by_block$`% over 75,000` + sj_nonworking_by_block$`percent less than 30` + sj_nonworking_by_block$`% speaking english > well` + sj_nonworking_by_block$`percent less than 1`)
summary(modeltest3)
## 
## Call:
## lm(formula = sj_nonworking_by_block$nonworking_percent ~ sj_nonworking_by_block$`% over 75,000` + 
##     sj_nonworking_by_block$`percent less than 30` + sj_nonworking_by_block$`% speaking english > well` + 
##     sj_nonworking_by_block$`percent less than 1`)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.9975  -3.2297   0.0295   3.1285  16.9204 
## 
## Coefficients:
##                                                    Estimate Std. Error t value
## (Intercept)                                        52.91517    3.21912  16.438
## sj_nonworking_by_block$`% over 75,000`             -0.10453    0.01486  -7.035
## sj_nonworking_by_block$`percent less than 30`      -0.07166    0.02927  -2.448
## sj_nonworking_by_block$`% speaking english > well` -0.04189    0.03161  -1.326
## sj_nonworking_by_block$`percent less than 1`       -0.07375    0.03210  -2.298
##                                                    Pr(>|t|)    
## (Intercept)                                         < 2e-16 ***
## sj_nonworking_by_block$`% over 75,000`             5.82e-12 ***
## sj_nonworking_by_block$`percent less than 30`        0.0147 *  
## sj_nonworking_by_block$`% speaking english > well`   0.1855    
## sj_nonworking_by_block$`percent less than 1`         0.0219 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.093 on 563 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:   0.19,  Adjusted R-squared:  0.1842 
## F-statistic: 33.01 on 4 and 563 DF,  p-value: < 2.2e-16

These variables do worse for the percent nonworking metric, which makes sense.